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Abstract 

We  extend  our  recent  work  on  higher-order  time  integration  of  Richards’ 
equation  to  layered  heterogeneous  porous  media,  using  a  differential-algebraic- 
equation-based  method  of  lines  (DAE/MOL)  approach.  We  show  that  the 
DAE/MOL  approach  is  robust  and  efficient  compared  to  standard  low-order 
time  integration  methods  for  heterogeneous  media.  We  also  show  the  ad¬ 
vantage  of  using  an  integral  representation  of  permeability  compared  to  a 
standard  arithmetic  mean  for  the  test  problems  considered  herein. 

1  Introduction 

Richards’  equation  (RE)  is  commonly  used  to  describe  flow  in  par¬ 
tially  saturated  porous  media  [] ,  although  questions  about  the  validity 
of  this  approach  remain  [].  RE  is  commonly  solved  using  low-order 
spatial  approximations  []  and  low-order  temporal-integration  meth¬ 
ods  Q;  several  codes  are  available  that  implement  these  methods  []. 

In  our  recent  work,  we  have  shown:  (1)  DAE/MOL  approaches 
for  solving  RE  are  robust  and  more  efficient  than  traditional  low- 
order  approaches  [?];  (2)  modification  to  a  standard  DAE  integrator 
can  further  improve  efficiency  [?];  (4)  use  of  vectorizable  interpolation 
methods  for  evaluating  constitutive  relations  can  drammatically  af¬ 
fect  simulator  performance  [?];  (3)  monitoring  of  the  condition  num¬ 
ber  of  the  Jacobian  is  an  effective  strategy  to  aid  in  the  selection 
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of  nonlinear  termination  criteria  [?];  (5)  lack  of  smoothness  in  van 
Genuchten-Mualem  relations  for  common  media  conditions  can  lead 
to  problems  with  convergence  of  nonlinear  solvers — a  situation  that 
can  be  remedied  using  an  integral  conductivity  representation  and 
cubic  spline  interpolation  [];  (6)  DAE/MOL  approaches  generalize 
to  two  dimension,  as  does  condition-number-based  termination  of 
the  nonlinear  solver  [];  and  (7)  transformation  approaches  can  sig¬ 
nificantly  improve  the  efficiency  of  RE  solution  methods  for  both 
standard  low-order  methods  and  DAE/MOL  approaches  []. 

While  these  advances  have  contributed  to  the  robustness  and 
efficiency  of  solutions  for  RE,  the  focus  of  this  work  has  been  on 
homogeneous  media,  except  our  recent  investigation  of  a  slightly  het¬ 
erogeneous  two-dimensional  system  in  which  the  focus  was  on  linear 
and  nonlinear  solver  issues  [].  Neither  our  previous  work,  nor  the 
work  of  others  that  we  are  aware  of  have  investigated  the  efficiency 
and  robustness  of  the  DAE/MOL  approach  for  solving  RE  for  het¬ 
erogeneous  media,  although  such  conditions  are  typical  of  natural 
systems.  This  seems  especially  important  to  consider  for  media  that 
lacks  smoothness  in  certain  situations. 

The  objective  of  this  work  is  to  compare  standard,  low-order 
methods  for  solving  RE  to  a  DAE/MOL  approach  for  heterogeneous 
media  conditions.  This  comparison  has  two  important  components: 
robustness,  or  the  reliability  of  the  solution  method;  and  efficiency, 
the  error  in  the  solution  achieved  for  a  given  investment  in  computa¬ 
tional  resources. 


2  Methods 

2.1  Formulation 


We  consider  two  common  forms  of  RE  in  one  spatial  dimension,  the 
pressure-head-based  form 


[c(V0  +  SsSa{^)]  ^ 
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and  the  mixed  form 


where  c  =  d9a/dip  is  the  specific  moisture  capacity,  Ss  is  the  specific 
storage  coefficient,  which  accounts  for  fluid  compressibility;  Sa  is  sat¬ 
uration  of  the  aqueous  phase;  ip  is  the  pressure  head;  t  is  time;  9a  is 
the  volumetric  fraction  of  the  aqueous  phase;  z  is  the  vertical  spatial 
dimension;  and  K  is  the  hydraulic  conductivity. 

We  consider  problems  with  auxiliary  conditions  of  the  form 
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where  Z  is  the  length  of  the  domain,  ipo  may  be  a  function  of  space, 
and  ip i  and  ip2  are  constants.  These  conditions  lead  to  the  devel¬ 
opment  of  a  sharp  infiltration  front  and  saturated  conditions  over  a 
portion  of  the  domain,  which  is  a  difficult  class  of  test  problem. 

For  closure,  we  use  the  standard  van  Genuchten  (VG)  pressure- 
saturation  relationship  [van  Genuchten(1980),  ],  which  is  given  by 


=  =  +  l^n,  ^<0 

Us  t/r  l  1,  Y  —  ^ 

and  the  Mualem  saturation-conductivity  relation  [Mualem(1976),  ] 

K  (Se)  =  KsS\l2  [1  -  (1  -  SlJn,-)mv]2  (7) 


where  rnv  —  1  —  1  /nv ,  Se  is  the  effective  saturation,  6r  is  the  resid¬ 
ual  volumetric  water  content,  9S  is  the  saturated  volumetric  water 
content,  av  is  a  parameter  related  to  the  mean  pore  size,  nv  is  a  pa¬ 
rameter  related  to  the  uniformity  of  the  pore-size  distribution,  and 
Ks  is  the  water-saturated  hydraulic  conductivity. 


2.2  Spatial  Discretization 

We  use  a  standard  finite-difference  approximation  to  discretize  RE 
with  respect  to  the  spatial  dimension  [Celia  et  al.(  1990) Celia,  Bouloutas,  and  Zarba, 
],  z,  where  z  E  [0,  Z].  We  consider  a  uniform  spatial  discretization 
comprised  of  nn  —  1  intervals  {[zj,^+ of  length  Az,  with 
Az  —  Zj{nn  —  1),  and  Zi  —  (i  —  l)Az  for  1  <  i  <  nn.  The  spatial 
operator 

MS") 
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is  approximated  at  z  —  z\  for  1  <  i  <  nn  by 

OsdiW  =  [Ki+1/2(tpi+1  -  tpi)  -  Ki_1/2{'ipi  -  4>i- 1)]  (9) 

+Az  1  [Ki+ 1/2  —  -ffj-1/2) 

where  nn  is  the  number  of  spatial  nodes  in  the  solution,  and  i/;,  is  the 
approximation  to  tp(zi). 

Interblock  conductivities,  are  evaluated  using  an  arith¬ 

metic  mean  (KAM) 


K-i±  1/2  —  {Ki  +  K%±l)  / 2 
and  an  integral  (KINT)  approximation 
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2.3  Temporal  Integration 

We  considered  two  temporal  integration  methods,  the  first  method 
is  a  standard  implicit  finite  difference  approximation  applied  to  the 
mixed- form  equation  and  solved  using  modified  Picard  iteration  (MPI), 
which  is  considered  the  standard  approach  [].  Time  step  size  was  se¬ 
lected  using  a  common  empirical  adaptive  time-step  control  algorithm 
[Rathfelder  and  Abriola(1994),  ] 

•  if  nn  <  mi  then  Afn+i  =  min(/tAf„,  Atmax) 

•  else  if  nn  >  mu  then  Afn+i  =  max(A tn/ft,  A tmin) 

where  m  is  the  number  of  iterations  required  by  the  nonlinear  solver 
to  converge  for  time  step  n,  mi  is  a  lower  iteration  limit,  rnu  is  an 
upper  iteration  limit,  ft  is  a  time-step  acceleration  factor,  A tmax  is 
the  maximum  allowable  time-step  size,  and  A tmm  is  the  minimum 
allowable  time-step  size. 

The  second  time  integration  method  investigated  was  a  DAE/ 
MOL  approach  for  the  pressure-head  form  of  RE  (1),  which  was  ap¬ 
plied  to  the  semi-discrete  form  of  RE  given  by 

=  °sdi ^ 
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where  A  includes  accumulation  and  compressibility  terms.  A  modi¬ 
fied  version  of  DASPK  [Brown  et  al. (1994)Brown,  Hindmarsh,  and  Petzold, 
Tocci  et  al.(1997)Tocci,  Kelley,  and  Miller,  ]  was  used  as  our  DAE 
solver.  DASPK  is  a  DAE  solver  based  upon  the  fixed  leading  coeffi¬ 
cient,  first-  through  fifth-order  backward  difference  formulas  and  con¬ 
taining  error  estimation  and  control  through  adjustment  of  method 
order  and  step  size  [Brenan  et  al.  (1996)Brenan,  Campbell,  and  Petzold, 


2.4  Efficiency 

We  define  efficiency  as  the  computational  effort  required  to  achieve  a 
specified  accuracy,  which  requires  evaluation  of  both  work  and  error — 
the  former  being  challenging  and  the  latter  trivial.  We  evaluate  work 
for  the  MPI  method  by 


Wp  —  wcnc  +  with 
and  for  the  DAE/MOL  method  by 

Wn  —  WjUj  +  wfrif  +  wini 


(13) 
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where  Wp  is  a  work  measure  for  MPI  methods,  wc  is  a  weighting  factor 
for  formation  of  the  coefficient  matrix  and  right  hand  side  vector, 
which  are  typically  done  at  the  same  time,  u>i  is  a  weighting  factor 
for  solution  of  the  linear  system  of  equations,  nc  is  the  number  of 
coefficient  matrix  formation  calls,  rq  is  the  number  of  linear  solutions 
performed,  Wn  is  a  work  measure  for  Newton  iteration  DAE  methods, 
wj  is  a  weighting  factor  for  formation  of  the  Jacobian  matrix,  wj  is 
a  weighting  factor  for  evaluation  of  the  function,  rij  is  the  number  of 
Jacobian  evaluations,  and  rif  is  the  number  of  function  evaluations. 
We  estimated  these  weights  based  upon  CPU  time  on  a  Hewlett- 
Packard  ...  running  version  ...  of  ...  operating  system,  and  using 
version  —  of  ...  complier  as  wc  —  0.484,  wj  —  0.552,  wj  —  0.271,  and 
wi  —  0.181.  Glenn  these  need  to  be  updated. 

Error  was  evaluated  by  comparison  to  dense-grid  solutions  using 
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where  k  is  the  norm  measure,  with  k  —  1  results  are  reported  in  this 
work;  and  ipi  is  an  approximation  of  the  true  solution  based  on  a 


Table  1:  Media  Properties 


Variable 

Sand 

Loam 

Clay 

9r  (  ) 

0.093 

0.078 

0.102 

0s  (-) 

0.301 

0.430 

0.368 

av  (m  x) 

5.470 

3.600 

3.350 

nv  {—) 

4.264 

1.560 

2.000 

Ks  (m/day) 

5.040 

0.250 

7.970 

Ss  (nr-1) 

1.0  x  10-6 

1.0  x  10-6 

0.000 

dense  spatial  grid.  The  dense-grid  solutions  were  generated  using  the 
DAE/MOL  approach  with  a  spatial  grid  size  equal  to  1/32  of  the  size 
used  in  the  test  simulations.  Glenn  check  this  too. 


3  Results 

Two  layered,  heterogeneous  test  problems  were  investigated  using 
the  MPI  and  DAE/MOL  approach  described  above.  For  both  test 
problems,  the  spatial  domain  was  z  E  [0, 4]  m.  The  temporal  domains 
were  t  E  [0, 0.25]  days  for  Problem  1  and  t  E  [0, 0.08]  days  for  Problem 
2.  The  media  for  Problem  1  consisted  of  four,  1-m  thick  alternating 
layers  of  loam  and  sand.  The  media  for  Problem  2  consisted  of  four 
1-m  thick  alternating  layers  of  clay  and  sand.  The  properties  of  these 
materials  are  listed  in  Table  1,  which  were  taken  from  the  literature 
[].  Problem  1  exhibits  a  larger  degree  of  heterogeneity  than  Problem 
2  due  to  the  wider  difference  in  Ks  and  nv  between  adjoining  lay¬ 
ers,  making  Problem  1  a  more  difficult  test  problem.  The  dense  grid 
solution  profile  for  both  test  problems  is  shown  on  Figure  1,  which 
illustrates  the  effects  of  media  heterogeneity  and  the  sharp-front  na¬ 
ture  of  the  solution  profiles. 

Quantitative  results  in  the  form  of  work  and  error  measures  for 
each  of  the  test  problems  are  shown  in  Figures  2  and  3.  The  compar¬ 
isons  of  the  MPI  and  DAE/MOL  approaches  showed: 

1.  the  DAE/MOL  approach  was  more  robust  than  the  MPI  ap¬ 
proach,  the  MPI  method  failing  to  produce  convergent  solutions 
for  Problem  1; 


Figure  1:  Solution  Profiles  for  Test  Problems  1  and  2. 

2.  the  DAE/MOL  approach  was  more  efficient  than  the  MPI  ap¬ 
proach  when  both  methods  converged,  producing  equivalent  ac¬ 
curacy  results  with  as  little  as  30%  of  the  computational  effort 
in  one  case;  and  Glenn  check  this. 

3.  the  KINT  produced  more  accurate  and  more  efficient  results 
than  the  KAM  approach.  Glenn  again  check  this. 
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Figure  2:  Solution  Profile  for  Test  Problem  1 
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Figure  3:  Solution  Profile  for  Test  Problem  2 
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